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r^ ; ABSTRACT 

O ,. We use weak lensing data from the Hubble Space Telescope COSMOS survey to measure the 

second- and third-moments of the cosmic shear field, estimated from about 450 000 galaxies 
with average redshift z ^ 1.3. 

We measure two- and three-point shear statistics using a tree-code, dividing the signal in 
jrt ' E, B and mixed components. We present a detection of the third-order moment of the aperture 

mass statistic and verify that the measurement is robust against systematic errors caused by 

point spread function (PSF) residuals and by the intrinsic alignments between galaxies. The 

amplitude of the measured three-point cosmic shear signal is in very good agreement with the 

predictions for a WMAP7 best-fit model, whereas the amplitudes of potential systematics are 

.,j. consistent with zero. 

Qv I We make use of three sets of large ACDM simulations to test the accuracy of the cosmo- 

T-j- ■ logical predictions and to estimate the influence of the cosmology-dependent covariance. 

• ' We perform a likelihood analysis using the measurement of {M^ ){9) and find that the 

^^ ■ flm — cTs degeneracy direction is well fitted by the relation: cr8(r2m/0.30)'''*^ — O.TS^q 2g 

f~^ ' which is in good agreement with the best fit relation obtained by using the measurement of 

We present the first measurement of the more generalised three-point shear statistic 
(A/f )(6'i, 92,03) and find a very good agreement with the WMAP7 best-fit cosmology. The 

X ■ cosmological interpretation of (Af3p)(6'i, 6*2, 6*3) gives cr8(f^m/0.30)°''^ = Q.69to°u- F™- 

S i thermore, the combined likelihood analysis of (Mf )(6'i, 6*2, 6*3) and {M^ ){9) improves the 

accuracy of the cosmological constraints to as{^m/0.30)'^'^^ — 0.69toi2j showing the high 
potential of this combination of measurements to infer cosmological constraints. 

Key words: Gravitational lensing - large-scale structure of the Universe - cosmological pa- 
rameters 



1 INTRODUCTION 

Weak lensing by large-scale structure, i.e. cosmic shear, constitutes a powerful tool to investigate the statistical properties of the dark mat- 
ter distribution. In the past years, the measurement of the observed correlation between pairs of galaxies has been successfully used to 
constrain the power spectrum of matter fluctuations (see Benjamin et al. 2007, Fu et al. 2008 and Schrabback et al. 2010 for the most 
recent results). With the amount of data available from current and future surveys, we should in principle be able to measure higher 
order statistics which provide a measurement of the non-Gaussian features of the matter distribution. In particular, third-order statis- 
tics are expected to increase the strength of the constraints obtained by weak lensing studies ( [Bernardeau, van Waerbeke & Mellier 1997[ 

* Based on observations made with the NASA/ESA Hubble Space Tele- operated by the Association of the Universities for Research in Astronomy, 

scope, obtained from the data archives at the Space Telescope European Inc., under NASA contract NAS 5-26555. 

Coordinating Facility and the Space Telescope Science Institute, which is t sembolon@strw.leidenunivnl 
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However, to date, predictions for three-point shear statistics are still inaccurate, as they are based on perturbative development and rely 
on the non-linear evolution of the power spectrum which is generally approximated by empirical models. Moreover, one needs to assess the 
accuracy of the measurement of the cosmic shear statistics in order to infer unbiased constraints on the cosmological parameters. In the past 
years, various publications focused on the study of two types of systematic effects which are likely to affect cosmic shear measurements. 
The first one is a measurement bias due to imperfections in the methods used to measure the ellipticity of galaxies and correct for the im- 
age point-spread-function (PSF) ( Heymans et al. 2006a , Massey et al. 2007a ; Bridle et al. 2009 , Kitching et al. 2009 , Semboloni et al. 2009J. 
The second one arises from the fact that galaxies align according to the local gravitational field ( Crittenden et al. 2001 .^Hirata & Seljak 2004| l 
creating intrinsic correlations. 

All these complications might explain why so far there have been only three attempts to measure three-point shear statistics: two using 
VIRMOS data (Bernarde au, van Waerbeke & Mellier 2003||Pen et al. 2003} and one using CTIO data ( |Jarvis, Bernstein & Jain 2004[ >. Due 
to issues mentioned above, neither of these papers presented a cosmological interpretation. 

Since then, there has been significant progress in the field and today we are able to quantify the impact of bias on the measurement 
of cosmic shear statistics, using simulations and various new diagnostics. Furthermore, the quality of the data has improved as well as our 
ability to interpret the measurement. 

In this paper, we present a measurement of three-point shear statistics obtained using the Hubble Space Telescope COSMOS data. We 
study the reliability of this measurement against systematics and use it to infer constraints on cosmological parameters. 

The COSMOS survey dScoville et al. 2007at combines high-resolution HST/ACS imaging covering a total area of 1.64 deg^ 
dScoville et al. 2007bt with deep ground-based multi-colour data, providing accurate photometric redshifts dllbert et al. 20091 1. Hence, this 
field provides a unique set of data to study the formation of structures and is particularly well suited for 3D weak lensing studies, such as the 
measurement of two-point shear statistics in tomographic bins (Massey et al. 2007b{ ^Schrabback et a l. 2010b . 

The paper is organised as follows. In Section |2] we introduce the notation and we test the accuracy of the cosmological predictions 
using various sets of TV-body simulations. In Section[3]we describe the data we use. In Section|4]we present the measurement of the two- 
and three-point shear statistics from COSMOS and test its reliability against systematic errors by PSF residuals and intrinsic ellipticity 
correlations. In Section[5]we explain the procedure used to determine the covariance matrix which we need to perform a likelihood analysis 
of our measurement. In Section|6l we perform a likelihood analysis and present the cosmological constraints we obtain. Finally, in Section|7] 
we discuss our results. 

2 THEORY 

2.1 Aperture mass statistics 

In this paper we present measurements of second- and third-order shear statistics obtained by using an aperture mass filter. Using the same 
notation as Schneider et al. (1998), the aperture mass is defined by: 

Map(6i) = j cFmeW^i'd) = f d^^QoW-fti^) , (1) 

where Ue is an axially symmetric filter of characteristic size 9 with zero mean. The second equality can be derived using the relation between 
the projected density field k and the shear field 'J — "ft + i^/x ■ We define the tangential 7t and the cross 7x shear components for a galaxy 
at the position ■»?, as the projections of 'y parallel and 45° rotated with respect to the line connecting the centre of the aperture and the centre 
of the galaxy. 

Aperture mass statistics are sensitive only to the E-modes of the observed shear field dCrittenden et al. 2002t . Replacing the tangential 
shear 7t in Equation (O with the cross component 7x one measures the non-gravitational, B-mode component of the observed shear field, 
Mx {9). A filter which is able to separate E- and B-modes unambiguously has a noticeable advantage when compared to other types of filters 
as this division can be used to test for the presence of systematics affecting the measurement of cosmic shear statistics. 

The relation between the two-point shear statistic {Allp){9) for a given characteristic scale 9 and the underlying cosmology is expressed 
bv dKaiser 1 9981 1 Schneider et al. I998t : 

{M!p){9) = 27r^n^ r" rfw4^ r sdsP(-^ w)[I{s9)f , (2) 

with 

g(w)= dwps(w) — , (3) 

where f-^lw) is the comoving angular distance, w is the radial coordinate, wh is the radial coordinate of the horizon; Ho, fim and ag are the 
Hubble constant, the matter density and the normalisation of the power spectrum of the matter fluctuations. (M^p) (9) depends on the power 

spectrum of matter fluctuations Pf , f . , w) and on the distribution of the sources ps(u;). Finally I{s9) is the Fourier transform of the filter 
A similar relation can be derived for (Affp) {9) in the perturbative regime dSchneider et al. 1998b : 

(Mfp)(e)= ° "g" / dw J \>. . / d;'s^p(-^,w)l{sr9) d''s2P(^^,w)l{s29)Ii\s-, + S2\9)F2{s,,S2) (4) 
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Figure 1. The top plots show the measurement of (A^f )(f?) for the set of simulations I (left panel): Om = 0.24, erg = 0.74 , II (middle panel): ilm = 0.25, 
as = 0.9 and III (right panel): Qm = 0.3, erg = 1. For set I we show the signal for two source redshift planes: z = 1.24 (black diamonds) and z = 0.84 (red 
triangles); for set II we show the signal for a source redshift distribution with (z) = 1.26 (black diamonds). For the the third set of simulations we show the 
signal for two source redshift planes: z = 1.12 (black diamonds) and z = 0.73 (red triangles). The three bottom plots show the measurement of (A/f ) (6) for 
the sets I (left panel), II (middle panel) and III (right panel) and for the same source redshifts. The theoretical predictions have been obtained using Equation 
(2) for {M^p){9) and Equation J4j for (A/fp}(9) for the same cosmology used to produce the simulations. The solid lines show the results obtained using 
the halofit model to compute the non-linear evolution of the power spectrum and the bispectrum. The dashed lines show the results using the non-linear power 
spectrum by PD96. The error-bars represent the dispersion between all the lines-of-sight available for each set of simulations. 



where F2(si, S2) represents the coupling between two density fluctuation modes characterised by the wave vectors Si and S2. 

One can estimate the value of Map (9) from the shear field as stated by Equation ([TJ, replacing the integral by a finite sum over the 
galaxies in the field. One can then use this measurement to estimate the second- and third-order moments of the aperture mass. However, on 
real data, this direct measurement does not allow one to separate the E- and B-modes unambiguously, as the properties of the filter Qgij}) are 
in principle lost due to the presence of masked regions. One can overcome this problem by estimating {M^p){6) and (A/fp)(S) using two- 



and three-point correlation functions as we will explain in more detail in Section MTT 



2.2 Testing the accuracy of the cosmological predictions 

In this section we investigate the accuracy with which {M^p)(6) can be predicted using Equation Q. Indeed, as we already pointed out, 
this equation is valid only in the perturbative regime so its precision is expected to decrease at small scales. We compute i<2(si, S2) using 
the fitting formula suggested by Scoccimarro & Couchman (2001). Scoccimarro & Couchman (2001) used N-body simulations from various 
ACDM cosmologies, and found that this fitting formula is able to predict the bispectrum within 15% accuracy. Van Waerbeke et al. (2001) 
used this same fitting formula to build third-order weak tensing predictions and found a discrepancy with simulations between 10% and 20%. 
The accuracy of the fitting formula is scale-dependent, thus we expect weak tensing predictions to have an accuracy which changes as 
a function of redshift, scale and filter type. For this reason, we directly compare {Mlp){9) predictions with values measured on three very 
large sets of ray-tracing simulations. All of them use a ACDM cosmology, and the main parameters of each set are summarised below: 

(i) The first set of simulations has been described by Vafaei et al. (2010); it is composed of 60 quasi-independent lines-of-sight, each 
having an area of 12.84 deg^. The cosmological model is very similar to the best-fit of the WMAP3 data l |Spergel et al. 2 007|>; JIa ~ 0.76, 
fim = 0.24, ilb — 0.04, h — 0.7 and as = 0.74. For each line-of-sight we have 40 convergence maps K{zi) corresponding each to a source 
redshift plane Zi. Each map contains 1024 x 1024 pixels of angular size 0.21 arcmin and all together the 40 maps cover a source redshift 
range between zero and three. 

(ii) As second set, we use the Millennium Simulation jSpringel et al. 2005]|Hilbert et al. 2009). The cosmology is characterised by JIa = 
0.75, flm = 0.25, fib = 0.045, ft = 0.73 and erg = 0.9. The set consists of 32 lines-of-sight of 16 deg^ each. Each map is 4096 x 4096 
pixels of angular size 0.058 arcmin. The redshift range is between zero and three. 

(iii) The third and last set of simulations, is described in Semboloni et al. (2007). The cosmological model is characterised by J^a — 0.7, 
flm = 0.3, h — 0.7 and erg = 1.0. The set includes 64 lines-of-sight of 49 deg^ each. Each map is 1024 x 1024 pixels of angular size 0.41 
arcmin. The redshift range is between zero and three. 
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Figure 2. Left panel: the red histogram shows the estimated redshift distribution of the faint sample. The green histogram shows the distribution of the 
individual photometric redshift estimates in the bright sample. The blue histogram shows the redshift distribution of the whole catalogue used to measure the 
two- and three-point shear statistics. The black line represents the best-fit to the total redshift distribution, obtained by using Equation (5) which is obtained 
for ZQ = 0.509, a = 1.29 and /3 = 0.95. Right panel: the black histogram shows the probability distribution of the best-fit value for zq in Equation (5) with 
a = 1.29 and /3 = 0.95. The probabihty distribution has been obtained producing one thousand bootstrap replications accounting for the estimated redshift 
error (see text). The blue line shows the Gaussian distiibution best-fit of the histogram. 

In general, we confiim that the fitting formula can reproduce the amplitude of {M^p){9) measured on simulations with the same degree of 
accuracy found by van Waerbeke et al. (2001). The goodness of the predictions depends on the angular scale (although for small scales we do 
expect the measurement to be affected by the finite size of the pixels) and slightly changes with the redshift. In the best cases the difference 
is about a few percent for {Mfp)(0) and ten percent for (Affp)(^). 

We show for example, in the top panels of Figure[T]the measurement of (Affp) (6*) using sets I (left panel), II (middle panel) and III (right 
panel) for various source redshifts. The error-bars represent the variance over the various lines-of-sight. The results have been compared with 
the predictions obtained using Equation ^ together with a non-linear power spectrum computed using the halofit prescription developed 
by Smith et al. (2003). Furthermore, we show models obtained using the fitting formula by Peacock & Dodds (1996) (hereafter PD96). 
The bottom panels show the measurement of {A/fp)(&) using sets I (left panel), II (middle panel) and III (right panel) for the same source 
redshifts. The results are to be compared with the predictions made using Equation ^ using either halofit (solid lines) or PD96 (dashed 
lines) to compute the non-linear power spectrum. 

In all cases (Affp) (6) predictions agrees fairly well with the measurements, although the models suggest values which are smaller than 
the measurements, especially for small angular scales. This is not the case for the set III and it can be explained by noticing that the pixel 
size for this set is too large to allow one to measure correlations at scales smaller than few arcminutes. 

For the two-point shear statistics the lack of agreement is due to the limited accuracy of the power spectrum in the non-linear regime. 
In this respect, our results are similar to the ones by Hilbert et al. (2009) who compared the projected power spectrum of matter fluctuations 
measured on the Millennium simulation with cosmological predictions by PD96 and by Smith et al. (2003), and found that both prescriptions 
strongly under-predict the amplitude of the power spectrum for small scales. 

We find that for {M^p){9) the disagreement is more significant especially at small scales. This is not surprising as the predictions depend 
on the non-linear power spectrum squared and they also rely on the perturbative approximation. The results obtained with the two methods 
are similar although PD96 underestimates the amplitude of the signal more than the halofit model. For this reason, we will use the halofit 
model throughout this paper. 

Overall, these results confirm an already known problem: in the near future we need to increase the accuracy of the cosmological 
predictions if we want to increase the accuracy of the cosmological interpretation of future data sets. This is not a problem for the aim of this 
paper which is meant to be a proof of concept and is still limited by statistical accuracy due to the small area of COSMOS. 



3 THE DATASET 

The HST/COSMOS Survey covers 1.64 deg^ consisting of 579 ACS tiles imaged in F814W with an exposure time of 2028s 
dScoville et al. 2007bt . A detailed description of our COSMOS shear catalogue, together with technical details of its derivation are given 
in Schrabback et al. (2010). Here we only briefly summarise its main features. 

Our galaxy shape measurements are based on the Erben et al. (2001) implementation of the KSB-l- formalism 
( [Kaiser, Squires & Broadhurst 1995[ [Luppino & Kaiser I997[ IHoekstra et al. 1998l l as adapted for space-based weak lensing measurements 
in Schrabback et al. (2007). We employ a principal component interpolation for the temporally and spatially varying ACS point-spread- 
function, parametric corrections for charge-transfer-inefficiency for both stars and galaxies, and a signal-to-noise dependent shear calibration 
correction. We match our shear catalogue to the ground-based photometric redshift catalogue from Ilbert et al. (2009), providing individual 
redshift estimates 2;poak and errors for 194976 galaxies with i^ < 25 (Subaru 3" Sextractor auto-magnitude). We will refer to this subsam- 
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Figure 3. Left panel: two-point correlations 5+ (9) (black solid line) and ^_ (6) (red dashed line) measured on the COSMOS galaxy catalogue. For l;^ (8) 
we show only the real component as the imaginary component is expected to vanish (^+(9) is real by construction). Right panel: amplitude of the three-point 
correlation functions F;, measured on quasi-equilateral triangles of side length 9. Because the bins are spaced in a logarithmic way triangles belonging to the 
bin s = t=|t — s|=6 are not strictly equilateral. The top-plot shows the amplitude of the real part, whereas the bottom plot shows the amplitude of the 
imaginary part. The en'or-bai's are computed as the dispersion measured in each bin and do not include sampling variance. 

pie as "bright sample". In addition, Schrabback et al. (2010) carefully estimated the redshift distribution pfaint (z) for 251 958 galaxies with 
1814 < 26.7 without individual redshifts, with a conservatively estimated ~ 10% uncertainty in the redshift calibration. We will refer to this 
subsample as "faint sample". In the left panel of Figure[2]we show the normalised redshift histograms for the faint and bright sample. 

In Section|6]we will present the likelihood analysis of the second- and third-order moments of the aperture mass measured from the total 
(faint-l-bright) sample. The amplitude of the signal depends on the redshift distribution of the sources Ps{z), thus the error on ps{z) affects the 
accuracy of the constraints. We show below that the error on the redshift distribution can be expressed by a single parameter; this parameter 
will be varied together with the cosmological parameters and it will act as nuisance parameter. 

We derive the redshift distribution of the total sample by adding the histogram of the bright sample and the probability distribution 
Pia.int{z) of the faint sample weighted by the number of galaxies contained in each catalogue. We fit the total redshift histogram using the 
following function JBrainerd et al. 1996l l: 



Ps{z) = 






^or(i±^) 



^"Pi-i;;; 



(5) 



In the left panel of Figure|2]we show the histogram of the total sample and the best-fit model, which is characterised by a = 1.29, /3 = 0.95, 
zo = 0.509 and average redshift z = 1.27. 

We now want to find the error on the parameters a, /3 and zo resulting from the individual redshift errors in the bright sample and 
uncertainty in the redshift extrapolation for the faint sample. We create one thousand bootstrap replications of the redshift distribution for 
the faint sample introducing a nuisance parameter fz similar to the one used by Schrabback et al. (2010). Using this parameter we create 
new redshift distributions piaintifzz), selecting fz from a Gaussian distribution centered at /z = 1 and having a dispersion a/^ — 0.10. 
We create one thousand bootstrap replications for the bright sample selecting the redshift value of each galaxy from a Gaussian distribution 
with dispersion az = {zf;s%,max — Z(;s%^rnin)/'2 where the quantities ^68%. min and 268%, maii define the redshift interval for which the 
probability distribution is higher than 68%. We build one thousand replications of the total catalogue by adding a bright sample histogram 
andapfaint(2) replicationof the faint sample. We notice that a and/3 are partially degenerate parameters of the Equation (O: fixing a = 1.29 
and f5 = 0.95 and then varying only zq allows for a degree of freedom which is high enough to obtain good fits for the various replications. 
Thus, we determine the best-fit parameter zq for each replication and build the probability distribution p{zo). We find that p{zo) follows 
fairly well a Gaussian distribution (see right panel of Figure[2} characterised by a standard deviation azg — 0.022, that is about 5% error on 
the mean redshift. 



4 MEASUREMENT AND ANALYSIS OF SYSTEMATICS 

4.1 Measurement of two- and three-point shear statistics 

The measurement of the two- and three-point correlation functions on the COSMOS field is accomplished by using a tree-code. The tree 
is built in a way analogous to the one described by Zhang & Pen (2005) and by Jarvis, Bernstein & Jain (2004). Each node of the tree is 
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Figure 4. The left panel shows the measurement of {M^ ){9) obtained using our COSMOS galaxy catalogue. The signal has been divided into gravitational 
(black diamonds) and non-gravitational (red triangles) component. The right panel shows the measurement of (M^ ){9) (black diamonds) and (MapM^ ) (9) 
(red triangles). The statistical noise affecting the E and B modes is plotted using solid-line eiTor-bars. The total noise on the E modes including sampling 
variance is shown using dashed-line error-bars. They have been computed using the set of simulations II . The measurement has been compared to a fiducial 
model obtained using a WMAP7 best-fit cosmology for the source redshift distribution of the shear catalogue (pink solid line). 



characterised by the average ellipticity (e) = (e^) + i{ey) over a region with centre xo, and characteristic size d. The characteristic size d 
is defined as the distance from the centre to the most distant galaxy included in the average. The step from a level to the next one is done by 
dividing the size d by two and recomputing the average quantities of the two new nodes; the new nodes are called "children". At the end of 
the division process we obtain a tree which has a top-level node ((e) , xo, d) computed over the whole image, whereas the bottom level has 
leaves (e, xo, 0) containing single galaxies; each node is linked to his children, the ensemble of children from the top node to a leaf is called 
branch. 

We compute the two-point correlation functions between &min ~ 5 arcsec and &max ~ 50 arcmin using logarithmic bins with size 
Slnd is 0.05. Starting from the top level, one must descend along the tree to compute the correlation functions between nodes of different 
branches. In order to decide whenever to continue the descent or to stop and compute the correlation functions between two nodes one must 
use a criterion. In our case we use the following criterion: ii±^ < S\n9 where di and d2 are the characteristic size of the nodes and s is 
the distance between the centres. When the condition is verified we compute the correlation functions between the two nodes and move to 
other branches. One can see that the smaller iSln^ is the deeper one needs to descent along the branches to verify the above criterion. When 
the criterion is satisfied we compute 



^+ = eie2*wiW2 and ^_ — eie2 exp 



•WiW2 



(6) 



where a is the angle between the x-axis and the line connecting the centres; ei, 62 are the average complex ellipticities and wi, W2 are 
the weights associated to each node, and • denotes the complex conjugate. We chose this weight to be equal to the number of the galaxies 
contained in the node. Averaging over all the pairs belonging to the same distance bin we obtain a measurement of the two-point correlations 
^+(e)ande-(e). 

For the three-point shear statistics we adopt a similar strategy, but now we have to consider that the correlations are function of three 
variables. Using the same notations as Jarvis, Bernstein & Jain (2004) we call s the vector connecting the closest pair characterised by sizes 
d2, da. The vector t connecting the second closest pair with sizes d2 and di, thus s < f < |t — s|. 

We compute the correlation functions for triangles with sides longer than Smin ~ 5 arcsec and smaller than Smax ~ 50 arcmin using 
again a logarithmic step with width Sln6 — 0.1. The criterion for the descent along the tree is similar to the one used for the two-point shear 
statistics 



di +d2 



< 5\n9. 



da 



± < Sine, ^1±^ < Sine 



It -si 



(7) 



When the above conditions are satisfied, we measure the four complex 'natural components' as suggested by Schneider & Lombardi (2003) 
and by Jarvis, Bernstein & Jain (2004). They are defined as follows: 



To = 616263 exp 



' U11U12W3 , Fi = 616263 exp 



' W1W2W3 , F2 = 616263 exp ''"' W1W2W3 , Fa = 616263 exp 



Wl'W2W3 



where a is the angle between s and the x-axis. 

We show in Figure [5] the amplitude of the two- and three-point correlations functions measured from the COSMOS galaxy catalogue. 
The left panel shows the amplitude of the correlation £,+ {e), and the real component of C-(^) as a function of the angular scale 6*; the 
top-right (bottom-right) panel shows the measurement of the real (imaginary) part of the four complex 'natural components' Fi measured on 
equilateral triangles as a function of the sides length e. 

In this paper we are mainly interested in three-point shear statistics and we could use the measurement of the correlation functions 
Fi directly to infer cosmological constraints ([Schneider, Kilbinger & Lombardi 2005). However, it is much more convenient to combine the 
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correlation functions to obtain the third-order moment of the aperture mass (Affp) {6) . In fact, not only the relation between (Affp) [6) and 
the bispectrum is much simpler than the one between the correlation functions and the bispectrum, but also, as we pointed out already, the 
aperture mass statistics allow one to separate the shear field in E- end B-modes providing a powerful way to detect systematics. 
By choosing 






0«('^) = TT^ e-P - ^ . Ue (^) = ^-^ (1 - 7-^) exp ^ — , (8) 



the relations between (Mfp)(6) and (A/fp)(6) and two- and three-point correlation functions are analytic dPen et al. 20031 ). With this partic- 
ular choice we obtain for the two-point shear statistics: 

{M^)(e) ^ {A4)(e)-{A/^)(e) + 2i{Af.pMx)(e) = |^e-Wr-(^), (9) 

{MM*){e) = (A/4)(e) + (M^)(e) = |^C+Wr+(^), (10) 

with 

In a pure cosmic shear field only the E-mode component (Affp) [9) is expected to be non zero whereas the B-mode component (Af ^ ) {0) can 
be used to investigate systematics, and the mixed component is expected to vanish if the non-gravitational modes are parity invariant. For the 
three-point shear statistics dPen et al. 2003l[Jarvis, Bernstein & Jain 2004| [Schneider, Kilbinger & Lombardi 2005[ ) the relations are: 

{M'){e)^{Ml^){e) + :ii{Ml^M^){0)-z{M^^Ml)(e)~i{Ml)(e) = &[ !^ f ^ro(s,t)ro(|,^) , (i2) 

{AfM*}{e) = (a4> {e) + ml^M^ ) {6) + {Kh^Ml ) {6) + i{Ml ) {6) 

M'Q f (ft r„ , s„ /s t\ „ , ,„ /s t\ ^ , ,„ /s t^ 

= '/ , , 

-^ s<t<|t — s| 

where 

rr u +\ q.f(if(if „ „ / g? + gj + g| \ ,,^. 

To(s,t) = — exp , (13) 



^L...,^ iih*''"^-(H)+r=i-"-"(i'l)+^''-')'-'&l)l 
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/(liflfcif gkaQa , qf + 2qfqf\ / q't + qi + qi' 



rr , ^.^ /qiq2 q3 i?iq2q3 , qi +2q2 qa \ / qi + <i2 + q-i\ ,,.. 

Ti{s,t) = -(^ — + jexp( ), (14) 



and 



qi = (s + t)/3s , q2 = (t - 2s)/3s , qa = (s - 2t)/3s , (15) 



and T2, Ta are obtained from Ti by a cyclic rotation of the indices. 

In the case of pure cosmic shear signal one expects to measure only the E-mode component (A4)(^)> whereas the two mixed com- 
ponents {M^pA4x)i6) and (AfapAfx)(^) and the B-mode component {Mx){0) are expected to vanish. If the observed field contains non- 
gravitational modes, the mixed component (Afap Af ^ ) (6) will not vanish whereas the two other components are still expected to vanish if the 
non-gravitational modes are parity invariant. Hence, (Afap A/x ) (^) can be used to quantify the amplitude of potential residual systematics in 
the galaxy catalogue. 

One notices that the filter in Equation {8]( has infinite support thus, in principle, one cannot compute (A/fp) (6') and {M'^p){9) unless the 
survey has infinite size. However, the filter has an exponential cutoff, thus in practice one can truncate the integration. Simulations show that 
to compute the second- and third-order moments of the aperture mass up to a scale 6 one needs to compute the correlation functions up to a 
scale ~ 46* 1 Jarvis, Bernstein & Jain 2004) ). We computed the correlation function up to scales of 50 arcmin and the second- and third-order 
moments of the aperture mass up to a scale of 12 arcmin. 

The left panel of Figure |4] shows the measured {Mlp){6) and {M^){6) after integration of the two-point correlation functions. The 
measured cosmological signal has been compared to a fiducial model obtained using a WMAP7 best-fit cosmology JKomatsu et al. 20T0t 
computed using Equation ^ with a redshift distribution given by Equation ^. The power spectrum in the non-linear regime is computed 
using Smith et al. (2003) and the transfer function introduced by Eisenstein & Hu (1998). 

The error-bars for {M^ ) (6) include only the noise due to the intrinsic shape of the galaxies. They have been derived using the noise maps 
computed as described in Section|5] We assumed that the error-bars on the non-gravitational modes are given only by shape-noise, although 
this might not be exactly the case as the B-modes might have a power spectrum if the PSF correction has some non-trivial residual pattern or if 
one accounts for intrinsic alignment effects. Assuming that the uncertainty is only the result of shape noise would underestimate the total error; 
nevertheless, as long as the amplitude of {M^){9) is much smaller than (A/fp)(S) the measurement is not significantly affected by system- 
atics. In our case, a linear fit performed using the shape-noise covariance matrix suggests a constant value {M^){9) = 0.13 ± 1.86 x 10~^, 
which is much smaller that the estimated cosmological shear signalj. For (Affp)(S) we show both the shape (solid error-bars) and the 

^ Note that Schrabback et al. (2010) used a different filter function for their from bootstrap resamples, leading to slight differences in both the signal and 

E-/B-mode decomposition using (Mf ) (6), for which they estimated errors errors compared to our results. 
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Figure 5. Left panel: the amplitude of the cosmological signal (Aff }(0) (black solid line) is compared with the measurement of (AfgYs)(^) ^^^'^ solid line) 
defined by Equation J17t which is expected to be non-zero for non-optimal PSF correction. For comparison we also show the amplitude of the aperture mass 
variance measured on the uncorrected stars, {M^g)(d) (blue solid Une). Right panel: the amplitude of the measured cosmological signal {AI^ ){9) (black 
soUd Une) is compared with the measurement of (M^ gYs ) (^) ^nd (Mj gYS )i^) (gr^Y ^nd red solid Unes) defined by Equations {TsJ and )19t . We plot for 
comparison the amplitude of the statistics (M^^^) [6) measured on the uncorrected stars. The error-bars on each statistic include only shape-noise. 



total noise (dashed error-bars). The right panel of Figure |4] shows the measurement of the third-order moment of the aperture mass com- 
ponents (A/fp) (61), {Ms,pM^){9) and the fiducial model computed for a WMAP7 best-fit cosmology using Equation l[4j and Equation ^. 
At small scales the amplitude of {M^p){9) is much larger than the one of the non-gravitational component {Mi^pM^){d). The amplitude 
of (AfapMx)(&) is consistent with zero, whereas the amplitude of the cosmic shear signal is in good agreement with the one suggested by 
the WMAP7 cosmology. At large scales the weak lensing signal becomes very small and thus its amplitude is comparable with the one of 
the mixed component, but still the amplitude of the mixed component is consistent with zero. The error-bars for {M^pMy){6) include only 
shape noise, whereas for {AI^p){9) we showed both shape (solid error-bars) and total noise (dashed error-bars). As we will discuss in more 



detail in Section|6]2] the sampling variance is the dominant source of noise for our set of data. 



4.2 Testing for systematics 

For a pure cosmic shear field one expects only E-type correlations between the observed ellipticity of galaxies. This means that whenever 
one finds B-modes it is likely that the observed alignment between galaxies is produced by something other than cosmic shear, for example 
PSF correction residuals. Therefore, the division of the two- and three-point shear statistics into E, B and mixed components shown in the 
previous section can be used to investigate the presence of systematic errors. 

However, one cannot expect the PSF residual to affect the E, B and mixed component in the same way. A more direct way to quantify 
directly the correlation between PSF and measured galaxy ellipticity is given by dBacon et al. 2003]|van Waerbeke et al. 2005) : 



esYs(e) = f^{e) 



(16) 



where ^ga is the correlation between the ellipticity of uncorrected stars, i.e. the PSF ellipticity, and the ellipticity of the PSF corrected galaxies. 
Here, ^ss is the correlation between the ellipticities of the uncorrected stars, ^sys(^) is normalised such that its amplitude can be directly 
compared with the amplitude of the cosmological signal, i.e. the correlation function ^{9) measured between pairs of galaxies. In fact, in the 
case of a perfect PSF correction the correlation between stars and galaxies would be zero so ^sys (9) would also be zero. On the other hand, 
if the PSF and the estimated shear are perfectly correlated the amplitude of ^sys(^) is about the same as (,{9). 

In a more general way, one can define cross-correlation two-point statistics for other types of filters such as the aperture mass. Integrating 
£,gs{9) and ^ss(9) one obtains the aperture masses {Mgg){9) and {Mss){9) and can define: 



(M|Ys)(e) 



(Mis) 



(e) 



(17) 



which is also normalised so that it can be compared with (A/fp) (9). 

For the three-point shear statistics we define the correlation functions: fgsal 
between two galaxies and a star, 5sss(^i, ^2, ^3) between three stars; here 6*1, 92, 



^1, ^2, ^3) between two stars and a galaxy, ^ggs{9i, S2, ^3) 
93 denote the distances between the vertices of the triangle. 
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Figure 6. Left panel: the amplitude of the cosmological signal {M^ ){9) expected for a survey with the same depth as COSMOS (black line) is compared 
with the amplitude of the expected |GI| (red line) for a survey with the same redshift distribution. The GI model has been computed using Equation j20t with 
the values Aqi = —1.29 X 10~^ h Mpc"'^ and 9qi = 0.93 arcmin which has been quoted by Heymans et al. (2006b) as the best fit-values for a model for a 
Universecomposedof about 30% elliptical and 70% spirals. Right panel: the amplitude of the expected cosmological signal (M^ )(9) (black line) is compared 
with Gil (red line) and |GGI| (green line). The values for Gil and GGI have been obtained by using the following values: Aau = 0.07 X lO"'^ h Mpc~^ 
and Sqqi = 3.74 arcmin, and Aqu = —0.22 X 10~^ h Mpc~^ and Oqqi = 1.61 arcmin, which are the best-fit models presented by Semboloni et al. 
(2008) for the same galaxy population used by Heymans et al. (2006b). 

Using these correlation functions we define the third-order aperture mass statistics (Mg^^ ) (9) , {M^^^) (6) and {Mf^^ ) (^) • Finally we define: 



{MlsYs){e) 



{MlsYs){0) 



{Ml)' 
(Ml.s) 



W, 



3/2 



(M|,)3/4 



■W: 



(18) 
(19) 



which are normalised so they can be directly compared with (Mfp)(S). 

The left panel of Figure[5] shows the amplitude of the cosmic shear signal, the amplitude of (Mfg)(S) and the amplitude of (71/|ys)(^)- 
The amplitude of {Mgg){6) is as high as the cosmic shear signal, whereas the amplitude of (M|ys)(^) is much smaller, meaning that the 
PSF is well corrected. The error-bars for {M^g){6) and (Af|Ys)(6) have been computed assuming Gaussian shape-noise both for galaxies 
and stars. The right panel of Figure [5] shows the amplitude of {M'ip){9), {M'iss){9), {A^i,sys)(^) and {-A^2,sys)(S)- It is very interesting 
to notice that the raw PSF does not show significant three-point correlation at least at small scales where the cosmological signal carries 
the most of the information. The PSF correction is very good at small scales thus, the cross correlation is always much smaller than the 
cosmological signal. 

One notices, that the definition of (A/j gYs)(^) and (Af|gYs)(^) is not unique. Indeed, one could have defined them by normalising 
{Mggs)iO) and {Mg^^){6) by a suitable power of {M^^^){9) instead than a power of {M^^){9). We find that the estimators (A/i^_sys)(^) 
and (A/|sYs)(^) defined in that way also suggest very low PSF residuals. Nevertheless, they are much more noisy than the ones shown in 
the plots of Figure|5]especially at small scales where the triplets of stars are very few. 



4.3 Intrinsic ellipticity correlations 

If galaxies are randomly oriented one can statistically interpret the observed correlation between galaxies as pure cosmic shear effect. 
However, galaxies belonging to clusters are affected by the local gravitational tidal forces. This effect creates an intrinsic alignment between 
galaxies which are physically close, so strictly speaking the hypothesis that galaxies are randomly oriented is not valid. Hirata & Seljak 
(2004) pointed out that the local tidal field affects the measure of two-point shear statistics through two terms: one of intrinsic alignments 
of physically associated galaxies and the other generally called shape-shear, which describes the correlation between the intrinsic shape of 
foreground galaxies with their local density field which is causing the shear of background galaxies. 

In the past years intrinsic alignment has been indicated as one of the potential sources of systematics affecting the measurement of cosmic 
shear statistics. Using TV-body simulations, Heymans et al. (2006b) found that for shallow surveys the intrinsic alignment can significantly 
affect the measured two-point shear statistics. This result agrees with the one that Mandelbaum et al. (2006) derived using galaxies from 
the SDSS spectroscopic sample. Semboloni et al. (2008) showed by using A'^-body simulations that the intrinsic alignment of galaxies can 
dominate the cosmological signal for three-point cosmic shear signal. However, according to these results the amplitude of the intrinsic 
alignment strongly depends on the projections and should not substantially bias the cosmic shear statistics measured on deep surveys such as 
COSMOS. 

The shape-shear term couples pairs of galaxies which are at different redshifts and it can in principle significantly affect the measurement 
of the two- and three-point cosmic shear statistics even for deep surveys. Heymans et al. (2006b) found that the correlation between the 
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cosmic shear ■y of a source at redshift Zs and the intrinsic shape e of a foreground galaxy at redshift zi close to the lensing overdensity can 
be modelled as: 

(7(..)e(.0)W = ^(^')/f---')/^';' (20) 

Jk(Zs) tl + Hgi 

where Aqi and Oqi are two parameters which depend on galaxy type. Semboloni et al. (2008) found using the same set of simulations that 
the shear-shape correlation affects the measurement of (Mfp)(0) through two terms that can be modelled by: 

GGI{zs,,Zs,,zue) = A(f^i - ^0 Mzs^ - ^0 ^^^^^i-o/Oao^) (21) 

Gll{z^,z,,9) = ^^lIlill^Giie'-"/*^") (22) 

fk(Zs) 

where the second equation uses the condition that two foreground galaxies that are close are both at the same redshift. Once again Aqgi, 
^Gii, ^GGi and 6*011 are parameters which depend upon the characteristics of the galaxies. Among the values of the parameters given by 
Heymans et al. (2006b) and Semboloni et al. (2008) we use the best-fit values as given for a mixed galaxy population to estimate the amplitude 
of the shear-shape coupling for the COSMOS survey; the mixed model is composed of about 70% spirals and 30% ellipticals and should 
be considered the most realistic among the ones analysed in those papers. Using the redshift distribution found in Section [3] we integrate 
Equations l l20b . l l21b . and ( 122b along the line-of-sight. 

The left panel of Figure |6] shows the amplitude of the expected GI component and the amplitude of the expected cosmological signal 
(Map) (6) for the same redshift distribution. The amplitude of the cosmic shear signal is always more than 20 times larger than the expected 
GI contribution. The right panel of Figure [6] shows the expected amplitude of the two components GGI and Gil and the expected cosmic 
shear signal {M'^p){9). As pointed out by Semboloni et al. (2008) the two components have different signs and the net effect depends on the 
redshift distribution. It turns out that for a deep survey such as COSMOS the contribution of the two terms is not only small but, the two 
terms almost cancel out so that the total effect from the shear-shape coupling is expected to have basically no impact on our measurement. 

Assuming that the empirical models we used to estimate the shear-shape are correct, the coupling is very likely not a source of significant 
systematic error for the measurement of the two- and three-point cosmic shear statistics. Nevertheless, we would like to point out that the 
effect of intrinsic alignment on two- and especially on three-point shear statistics is still poorly known and generally based on approximations 
and toy-models which may be too simplistic to describe reality. Much more work is needed to rigorously quantify the effect of intrinsic 
alignment on cosmic shear statistics. 

5 COVARIANCE MATRIX ESTIMATION 

The covariance matrix C can be naturally divided into three components JSchneider et al. 2002b : 

where Css is the pure sampling variance term, Cnn is the statistical noise, which we assume is generated only by the intrinsic ellipticity of 
the sources, Csn which is a cross term generated by the correlation between cosmological signal and statistical noise. 

To estimate the covariance matrix for the two-point shear statistics one can use analytical approximations dSchneider et al. 20021 
IJoachimi, Schneider & Eifler 2008^ . However, analytical approximations assume that the field of matter fluctuations follows a Gaussian 
statistic; this causes the covariance matrix to be underestimated at small angular scales iSemboloiri et al. 2007al ITakada & Jain 20091 
IPielorz et al. 20 09). The third-order shear moments are generated by the non-linear evolution of the density matter field, i.e. when the 
Gaussian approximation is not valid, so it is even less desirable to use the Gaussian approximation to compute the associated covariance 
matrix. For this reason we compute the covariance matrix using ray-tracing from the Millennium Simulation (i.e. the set simulations II). 
We divide each of the 32 fields into 9 subfields having the same area as COSMOS, thus obtaining 288 quasi-independent lines-of-sight. 
For each subfield we generate shear catalogues with the same number of galaxies, the same redshift distribution and the same masks as the 
COSMOS catalogue; finally we need to add shape-noise to each simulated galaxy. To do so we need to determine the ellipticity dispersion 



j^ + o-g 2 of our catalogue. In order to investigate if we need to vary a^ with redshift, we divide the catalogue in magnitude and 
size bins and inspect the variance in each bin. Note that the eventual change of tJe is not necessarily the result of galaxy evolution, but it is 
mainly due to the fact that the shear estimation is more affected by noise for faint and small objects. 

We find that the dispersion depends on the magnitude; for bright objects the dispersion is a^ = 0.39, for faint objects ae — 0.49. 
However, this is only a minor effect and should barely change the value of ae as a function of the redshift (the very bright galaxies have a 
dispersion 20% smaller than the faint ones). Thus, we decide to use a dispersion ae = 0.44 average of the whole sample, overestimating 
the noise at low redshifts. This should not have a significant effect as we expect the sampling variance to be the main source of the noise. 
Furthermore, we know that most of the contribution to the weak lensing signal comes from sources at high redshift. Once we added the noise 
to the shear catalogues, we compute the two- and three-point correlation functions on the 288 lines-of-sight and integrate to obtain (M^p) (9) 
and(A'4)(e). 

Finally, we compute the covariance between two measurements Xi and Xj (where Xi can be either the measurement of {Affp)(S) or 
(Mfp) (9) at a given angular scale) as: 



*~'\Xi^Xjj — y yXi^n '^^}\^j 



(24) 
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Figure 7. Left panel: probability distribution for the parameters Vim and cs obtained preforming a likelihood analysis of the measurement of (M^ )((?). 
The two different colours represents the regions where the confidence level is higher than 68% and 95%. Middle panel: same as left panel obtained using 
the measurement of {M^ ){8). Right panel: probability distribution for the parameters Qrn and trg resulting by a likelihood analysis in which we used the 
measurement of (Mfp)(e) and (A/fpXe) together The solid (dashed) lines represent the 68% (95%) level of confidence for (Affp>(e) and (J\/|p)(e) used 
separately as in the left and middle panels. 

where Xi is the measurement of the quantity Xi averaged over the 288 lines-of-sight. 



6 COSMOLOGICAL CONSTRAINTS 

6.1 Likelihood analysis 

The measure of two- and three-point shear statistics can be used to infer cosmological constraints on ACDM models. This is done through a 
likehhood analysis which establishes the posterior probability distribution P(pjx) for the ensemble of parameters p: 



P(p|x)(x£(x|p)P(p), 



(25) 



where P(p) is the prior probability distribution for the ensemble of parameters and x is the vector of n measurements. We will assume that 
the likelihood £(p|x) is Gaussian: 



£(x|p) oc exp 



¥- 



x™.(p))*C ^(x-x,„(p)) 



(26) 



although this is only an approximation, whose limits have been discussed in Hartlap et al. (2009). The vector x contains the measurements 
of {Mfp)(6l) and {M'ip){9) presented in SectionE] Xm(0) is the vector of the models of {Ml^){e) built using Equation l|2} and {Ml^){9) 
built using Equation Q; C~^ is the inverse of the covariance matrix computed in Section|5] Before inverting it we multiply by the correction 
factor: 



{r. 



in„ 



(27) 



which allows one to obtain an unbiased estimation of the inverse of the covariance matrix for a given number of simulations rimaps and a 
number of measurements n ( [Hartlap, Simon & Schneider 2007| l. 

As the COSMOS field covers a small area we do not expect to have strong constraints if we allow a large number of cosmological 
parameters to vary. Hence, we decide to vary only the matter density Qm and the power spectrum normalisation erg to which the cosmic shear 
statistics are particularly sensitive. We assume uniform priors for both parameters; we vary Q,,n between [0.1, 1.] and erg between [0.5, 1.2]. 
Moreover, we found in Section[3]that the redshift uncertainties can be included in the likelihood analysis by marginalising over the parameter 
zo between [0.443, 0.574], keeping the other two redshift parameters fixed: a — 1.29 and /3 — 0.95. This interval in zq represents the ±3cr 
width of the Gaussian distribution of p(Zo) found in Section[3] We fix all the other cosmological parameters to the best-fit cosmology of 
WMAP7. 

The left and middle panel of Figure [7] show the likelihood contours in the [Qm, erg] space, marginalised over the zo parameter, obtained 
using {Mlp){6) and (Affp) {6) respectively; the right panel shows the likelihood contours obtained by combining the measurements of both 
statistics. For each of the posterior probability distributions in the [Sim , fg] space after marginalisation over the redshift parameter zq we find 

the local maxima and fit them using a relation erg ( jj-gj J = Y where X and Y are free parameters. This allows us to fit the degeneracy 
direction between Q.m the ag fairly well obtaining the following results: 



a"8 



0-8 



0.30 y 



*1 0.30 



0.70irii (M^p)(e) 



= 0.78t°;^^ {M^p)(e) 



= 0.69t«;li (M^p)(e)and(M^p)(0) 



(28) 
(29) 
(30) 



where the la error-bars correspond to the upper and lower limit of the 68% confidence regions. 
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Figures. The upper row of plots shows the logarithm of the amplitude of the (Mf )(6) covariance matrix for the six angular bins we use for our cosmological 
interpretation. From the left to the right we show the results for the set of simulations I, II and III. The bottom row shows the logarithm of the amplitude of the 
{M'^ ){0) covariance matrix. In all cases we rescaled the amplitude of the covariance matrix to the size of the COSMOS survey as indicated in the text. 

As expected, the constraints on the cosmological parameters, especially the ones obtained using {M'^p){0) are broad as the COSMOS 
field is small, and the cosmic variance is large. Moreover, the measurement of {M^p){9) does not improve the accuracy of the cosmological 
constraints when used together with the measurement of (Affp) (6). This is not surprising, as Vafaei et al. (2010) pointed out that the geometry 
of a narrow and deep field such as COSMOS is not optimal for the measurement of three-point shear statistics. According to their results 
wide and moderately shallow surveys would be more suitable and would allow one to partially break the degeneracy between Qrn and as by 
using two- and three-point shear statistics. 

We notice that our result for (Map) (0) is in agreement at the la level with the result ag ( q-^ ) — 0.81 ± 0.17 found by Massey et al. 
(2007b) who previously measured two-point shear statistics on the COSMOS catalogue performing an analysis which is completely indepen- 
dent from the one presented in this paper Our results are also in excellent agreement with the 2D analysis result as [ g-gj 1 = 0.68 ± 0.11 
presented by Schrabback et al. (2010) who used the same catalogue. The very small difference between the results may be originated from 
the fact that Schrabback et al. (2010) used the correlation function to perform a likelihood analysis and a different type of binning which 
effectively probes the power spectrum at slightly different scales. In addition, they quoted marginalised means, whereas our fit describes the 
maximum posterior, which typically corresponds to slightly higher as ■ 

The likelihood analysis of (Mfp) {9) gives contours which are fully consistent with the ones obtained using {M^p){9), although they are 
broader. The best-fit cosmology for {M^p){9) and {Mlp){9) is similar suggesting that our measurement is not affected by large systematics 
which, in principle should affect the two- and three-point statistics in a different way. 

To investigate the effect of poor modelling on the cosmological constraints we just derived, we repeat the likelihood analysis using 
as data vector the measurement of {M^p){9) and (Map)(^) measured on the Millennium ray-tracing simulations (i.e. set II). For this test 
we change our WMAP7 priors to the ones in the Millennium simulation. We find that for {M^p){9) the best-fit is as = 0.92lgj2 for 
n,„ = 0.25, whereas for {Mlp){9) the best-fit is erg — O.OS^Qjy. The difference between the best-fit values of erg obtained from the 
likelihood analysis of {M^p){9) and {Mlp){9) is marginal and the inaccuracy of the cosmological predictions should not bias significantly 
our results. 

6.2 Effect of the cosmology on the likelihood analysis 

The result in the previous section have been obtained using the covariance matrix from the Millennium Simulation. It has been pointed out by 
Eifler, Schneider & Hartlap (2009) that in the case of two-point shear statistics, the cosmology used to compute the covariance affects the error 
on the cosmological constraints. In fact, we expect the effect to be even stronger for the three-point shear statistics given their sensibility to 
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Figure 9. Left panel: noise-to-signal ratio expected for the (Aff )(6) measurement on the COSMOS survey, estimated using the set of simulations I (solid 
lines) II (dashed lines) and III (dot-dashed lines). The solid black lines show the total noise-to-signal ratio, whereas the green lines show the amplitude of the 
noise-to-signal coming from sampling variance, the red lines show the noise-to-signal given by shape noise and the blue lines show the noise-to-signal from 
the cross term. Right panel: same as left panel for the measurement of (A/f )(6'). 



(jg. Eifler, Schneider & Hartlap (2009) suggest to compute the likelihood in each point of the parameter space by using the covariance matrix 
obtained from the cosmology characterised by those same parameters. This would require us to have large sets of ray-tracing simulations 
for each point in the parameter space we want to explore, which we do not have. However, since we have three sets of simulations which 
explore a quite large range of erg values, we can use all of them to compute covariance matrices and evaluate the differences in the likelihood 
analysis. For this test, we compute the covariance matrices using the direct measurement of the aperture mass statistics. The covariance 
matrices obtained in this way reflect less accurately the noise specific to our survey and our estimator, but they can be computed much faster. 
For each set, we create lines-of-sight weighting the convergence maps K[zi) according to Pb{z) given by Equation l|5j to have the same 
source redshift distribution as the best-fit to the COSMOS data. We account for intrinsic ellipticity dispersion by adding to each pixel of the 
corresponding convergence noise defined as dvan Waerbeke 2000llVafaei et al. 2010) : 

^2 

(31) 



2n6»2 



where iTe = 0.44, 6 is the pixel size and n = 75/arcmin'^ is the density of galaxies. 

From Equation ^ one notices that the filter C/g(i9) has infinite support thus the measurement of Ma_p{6) on a survey of finite size 
is biased. In practice, for a given characteristic size 9 we can truncate the measurement of Afap(S) at a scale 49. Because of that, given a 
simulation of side length I the effective area used to compute the aperture mass is (l — 86)^ as one cannot compute the aperture mass in 
regions which are closer than 46 to the borders. Hence, the effective area used for each bin 6 is different, leading to a covariance matrix 
which has noise properties which are different for each bin and it depends in a complicated fashion on the angular scale and on the area of the 
simulations. In order to avoid this issue we do not compute the aperture mass in regions which are closer than 49max to the borders, where 
dmax is the maximum angular bin we measure {Mlp){9) and {Mlp){9). With this cut the effective area used for the simulation is 4.12 deg^ 
for set I, 5.99 deg^ for the set II, 29.68 deg^ for set III, whereas the initial area of each set was: 12.84 deg^, 16 deg^, 49 deg'^. As one can 
see the cut we applied reduces the area of each set of simulation considerably. This increases the noise affecting the covariance matrix but 
allows us to have an homogeneous measurement of the covariance at all scales using (almost) the same area of the sky. 

Once we have the final noisy convergence maps, obtained simply by adding the shape-noise maps to the convergence maps we estimate 
for each pixel the value of the aperture mass using the relation between the convergence k and the aperture mass Afap in Equation ([TJ. Finally 
we compute for each field the second- and third-order moments. We can rescale the final covariances multiplying by the ratio of the effective 
area used for each set of simulations and the area of the COSMOS survey (Schneider et al. 2002). 

In Figure[8]we show the amplitude of the elements of the final covariance matrix both for (A/fp)(6) and (A/fp)(S) in the six angular 
bins we use for the cosmological interpretation. As one can see the overall shape is the same for the three simulations. However, as one 
would expect, the overall amplitude of the covariance matrix is larger when the value of erg is higher. In Figure|9]we show the noise-to-signal 
ratio for each of the terms in the covariance matrix we just computed, i.e. the amplitude of the diagonal part of ^/C7I, \/Cs„ and \fC^ over 
the cosmic shear signal for the set of simulations I, II and III. The plots indicate that for a narrow and deep survey as COSMOS the error 
generated by sampling variance Cas is always dominant, as the contribution of statistical noise is minimised by the high density of galaxies. 
In particular, the amplitude of the diagonal part of \fCTs for {M^p){9) is already as large as the signal at small scales. 

Figure [To] shows the size of the likelihood contours obtained for each set of simulations for {Allp)(6) (left panel) and {M^p){9) (right 
panel). For each likelihood analysis we used as a data-vector a fiducial model with the same cosmology used to create the simulations. One 
can see that especially in the case of (Affp) (6*) the contours are similar. The signal-to-noise does not depend strongly on the cosmology. 
This implies that the covariance matrix strongly depends on the cosmology. Indeed, the sampling variance term depends very roughly on 
(jgSlm as it is proportional to the signal squared (compare Equation|4l(. Furthermore, a higher erg value indicates that, for the same redshift 
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Table 1. The table on the left shows the best-fit values and ±1(t error- 
bars for the erg parameter for a fixed Qm =0.3, obtained performing 
the same likelihood analysis for the COSMOS data but with covariance 
matrices derived using different cosmologies. For this comparison we 
used covariance matrices which have been computed directly from the 
moments of the aperture mass statistics (see Section l6!2t . Notice that 
because of that the ±1(T error-bars on ag are slightly different from the 
ones obtained using the covariance matrix computed from the correla- 
tion functions. 
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Figure 10. The left plot shows the one sigma constraints in the [0,n,crg] space obtained for a 20 deg'^ simulated survey with the same depth as COSMOS using 
(Aff ) (0) and the covariance matrix from the set of simulations I (blue solid contours), II (red dotted contours) and III (pink dashed contours). To perform the 
three likelihood analyses we used the fiducial model which corresponds to the cosmological model employed to produce each of the set of simulations. The 
right plot shows the same results for {M^ )(9). 

distribution, the structures are more evolved thus the modes are more coupled and the measurements are more correlated. Since the choice 
of the fiducial cosmology has a large impact on the cosmological constraints, therefore one should choose a fiducial cosmology as realistic 
as possible. We quantify the impact of the choice of the covariance matrix in Table showing the error-bars one would obtain performing a 
likelihood analysis of the COSMOS data using the covariances I, II, III. A posteriori, one can see that it is reasonable to use the set I or II to 
interpret our data as the observations rule out the cosmology of the set III. Furthermore, our analysis confirms that it is important for future 
missions, such as Euclid and IDEM, to be able to perform a likelihood analysis, varying the covariance matrix according to the parameters 
as already suggested by Eifler, Schneider & Hartlap (2009). 



6.3 Measurement of (A/fp)(6ii, 6I2, 6I3) 

In the previous section we presented the cosmological constraints obtained measuring (Affp) (S) on the COSMOS field and we commented 
that these constraints are weak due to the small size of the data sample. However, we ignored that the measurement of {M^p){6) contains 
only one part of the information about the shape of the bispectrum. Indeed, as it has been pointed out by Schneider et al. (2005), the amplitude 
of {Ailp){6) mostly depends on the shape of the bispectrum B(|ki|, |k2|, |k2 — ki|) for modes with |ki| ~ jk2| ~ |k2 — kij ~ 1/6, 
i.e. on equilateral triangles. One can easily generalise Equation Q using filters with different characteristic size defining (M^p){di, 02,03). 
By varying independently the three characteristic filter scales one explores the amplitude of the bispectrum measured on a larger variety 
of triangles. In this section we aim to improve the accuracy of the cosmological constraints obtained in the previous section measuring 
{M^p){0\, 02,03) varying 0\, 62 and 63 independently. Each angle can take six values between 0.8 arcmin and 12 arcmin with the condition 
Oi ^ O2 ^ O3 as the measurement is symmetric toward any permutation of the filter scales. In this way we obtain 56 measurements of 
(Mfp) {61 ,02,03). Table 2 presents the way in which we order the various triplets of filter scales and the left panel of Figure [TT] shows the 
signal we measured for the triplets. In order to check the cosmological origin of the signal, we divided the signal into E-, B- and mixed-modes. 

The amplitude of the E-modes signal is in very good agreement with the predictions. Moreover, the amplitude of the non-gravitational 
component {MapAfx)(^i, 02,63) is generally low and consistent with zero in most of the bins, confirming that this estimator is not substan- 
tially affected by residual PSF systematics and can be used to constrain cosmology. 

Ideally we would like to include the whole set of new measurements in the likelihood analysis. However, we have a small set of 
simulations, so the covariance matrix may be poorly conditioned and lead to biased constraints if we choose too many measurements or 
measurements which are highly correlated. 

The 56 X 56 covariance matrix has a condition number, i.e. the ratio between the highest and the smallest eigenvalues, of 10^^ indicating 
that as long as we keep all the scales the inversion is unstable. In fact, with such a high condition number every modelling error will have a 
huge impact on the likelihood result. Note that since we use a logarithmic step, the 56 triplets are mainly composed by triplets where Oi, 62 
and 63 are similar and small. These triplets explore the bispectrum at very similar modes. Hence, their measurement is highly correlated. 
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Table 2. The first three columns show the values taken by the 
first, second and third aperture mass angles: 8i, 82, 63, where 
[9i, ■ ■ ■ , 6e] indicate six angular logarithmic bins between 0.8 
arcmin and 12 arcmin. The last column shows the number as- 
sociated to each triplet shown in the left panel of Figure [TT] 



In order to reduce the number of triplets having similar filter scales, we retain only triplets with 61 — 62 if 61 < 4 arcmin. If 61 > 
4 arcmin the filter sizes are very different so the triplets can be all retained. With these conditions we are left with 25 triplets. The condition 
number for the new 25 x 25 covariance matrix is 10^, indicating that the likelihood analysis with these triplets is well-conditioned. Hence, 
we perform a likelihood analysis using these 25 {M^p){9i,92, 63) measurements from the COSMOS data. We show in the middle panel 
of Figure[TT]the posterior probability distribution for Q,n and as after marginalisation over the redshift parameter zq. One can see that the 
cosmological constraints are significantly improved. The best-fit as — 0.69'^q'^I^ for Qm ~ 0.3 is in good agreement with the constraints we 
found using (Mfp) (0) and (A/fp) (6) . In the right panel of FigurefTTI we show that the combined analysis of (Mfp (61) and (Mfp) {9i ,62,03) 
increases the strength of the cosmological constraints. Indeed, the combined likelihood results yields as(Qm/0.30) '^ = 0.69to'i2- 

Using Fisher matrix analysis, various authors ( [Kilbinger & Schneider 2005] ITakada & Jain 2005] |Berge, Amara & Refregier 20I0[ | 
showed that (A/fp)(Si, 62, 63) contains more information than {M^p){9). They also suggested that since {M^p){6) and {M^p){9i, 62,03) 
depend differently on Q,m and erg they can be combined to partially break the erg — 0,m degeneracy. In this respect our results are in good 
agreement with their prediction and show the potential associated to this measurement. 



7 CONCLUSIONS 

In the past decade several papers discussed the use of the three-point statistics to explore the non-linear features of the field of matter 
fluctuations. It has also been pointed out that adding three-point to the two-point shear measurements allows one to increase the accuracy 
of the constraints on cosmological parameters (see Vafaei et al. 2010 for a quantitative analysis based on recent ray-tracing simulations). 
However, measuring three-point shear statistics appears much more challenging and only few attempts have been done so far. The aim of 
this paper was to show that today we have the expertise and the data quality to allow for the measurement of third-order cosmic shear 
statistics. For such a purpose, we used the COSMOS data set: a mosaic of contiguous pointings of space-based data which represents the 
ideal benchmark to test the quality of weak lensing analyses for future space-based data. In order to show the reliability of our results we 
focused our attention on some very important aspects to correctly interpret three-point shear statistics. 

As first step, we tested the accuracy of the cosmological predictions against four sets of recent ray-tracing simulations based on ACDM 
cosmologies. We found that current cosmological predictions for the two- and three-point shear statistics used in this paper, i.e. {M^p){6) 
and {Mlp){0), are able to reproduce the signal measured in the simulations typically only within 10% and 20% accuracy respectively. This 
clearly demonstrates that improved model predictions will be required for the analysis of future surveys. 

We used empirical models derived using A'^-body simulations to quantify the effect of the shape-shear coupling on the measurement of 
the two- and three-point shear statistics. We found that for the COSMOS Survey this effect should be subdominant. To date, the shear-shape 
effect is still poorly constrained and the empirical models we used are based on a toy-model which might not be realistic enough. More work 
still needs to be done to understand and model the effect of the intrinsic orientation of galaxies on cosmic shear statistics, although the recent 
results by Mandelbaum et al. (2009) suggest that the shear-shape should affect the value of erg from two-point shear statistics at most at the 
2% level. 

In order to demonstrate the robustness of our measurement, we divided the signal into gravitational and non-gravitational components 
and investigated the amplitude of the cross-correlations between galaxies and stars, which can be used to reveal the existence of PSF residual 
correlations. We found that the measured (A/fp)(S) signal has an amplitude which is in very good agreement with the WMAP7 best-fit 
model, and that all conducted tests are consistent with negligible residual PSF systematics. Moreover, we showed that the non-gravitational 
component is consistent with zero both for the two- and the three-point shear statistics. The cross-correlation between stars and galaxies is 
also consistent with zero. 

The likelihood analysis of (A/fp)(6) gives as = O.TSljJ'jg for fixed flrn ~ 0.30. The likelihood results obtained using (A/fp)(S) alone 
are in very good agreement with the ones obtained using {Ai'lp){6). Our error budget is dominated by the large sampling variance which 
strongly affects the measurement of the three-point shear statistics on a small field such as COSMOS. In fact, joining the two measurements 
we obtain the same cosmological constraints derived by using {Mlp){6) alone. This is in agreement with the results by Vafaei et al. (2010), 
who found that a combined measurement of {M^p) (6) and {M^p){6) on narrow deep surveys such COSMOS cannot be used effectively to 
break the [Sim , f-g] degeneracy. Finally, we generalised the definition of ( A/fp ) (^) to ( A/^p ){di, 62,63) and we showed that its amplitude is in 
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Figure 11. Left panel: the black diamonds show the amplitude of the cosmological signal (A/f )(9i, 82,03) measured for each of the 56 combinations of 
61, 82 and 8-j as a function of an identification number which we assigned to each triplet (see table[2)- The amplitude of the cosmological signal is compared 
with the fiducial WMAP7 cosmology (pink solid line) and the non-gravitational component (A/apM^ ){9i , 82, 83) is also shown (red diamonds). The solid 
error-bars for (Mf ) {81, 82, 83) include only shape-noise and we omitted the identical error-bars for the non-gravitational component for clarity. The dashed 
error-bars represent the total errors computed using the ray-tracing set II. Middle panel: likelihood analysis contours obtained using the measurement of 
(-^ap)(^li^2, 93). The black solid (dashed) contours show the likelihood 68% (95%) contours obtained using {M^p){8) for comparison. Right panel: 
likelihood analysis contours obtained using the measurement of (Mf )(0i, 82,83) combined with (M^ )(8). 



excellent agreement with the WMAP7 best-fit model. We find, as expected, that this measurement improves the accuracy of the cosmological 
constraints ag = 0.691q54 for fixed Q,m = 0.30. The combined analysis of {M^p{9) and {M^p){6i, 62,03) further increases the strength 



of the cosmological constraints. Indeed, the combined likelihood results yields as{ilm/0.30)'^'^'^ = O.GOIq'jI- 

We find this result very encouraging, however we would like to make clear once more that the lack of precision which affects the non- 
linear evolution modelling, is today still significant. We already discussed how the cosmological constraints we obtained are sensitive to the 
non-linear evolution of the power spectrum. Note that Lawrence et al. (2010) very recently provided a substantially improved prescription for 
the non-linear power spectrum , but similar results will also be required for the bispectrum. Various aspects which have been neglected in this 
paper are likely to affect our results. For example, we did not account for massive neutrinos and the effect of baryons on small angular scales: 
both are expected to change the evolution of the matter fluctuations in particular in the non-linear regime. Furthermore, we did not account 
for lens-clustering effects which are also expected to change the amplitude of the measured three-point shear statistics dHamana et al. 20021 . 
Constraints from wide surveys will more strongly benefit from the inclusion of third-order statistics. Hence, for future surveys such 
as Euclid, LSST and IDEM covering thousands of square degrees it will be possible to use the measurement of the two- and three-point 
shear statistics together with photometric redshift information to infer tight cosmological constraints. In this paper we did not make use of 
the photometric redshifts estimated for individual sources; we expect the tomographic measurement of the third-order shear statistics to add 
further information and improve the precision of the cosmological constraints. 

In the prospective of these future weak-lensing missions, it is important to show that it is today possible to measure three-point shear 
statistics, to quantify the level of systematics and to interpret the cosmological origin of the signal. It is also important to point out which 
are the limiting factors requiring further work. This paper shows that the measurement and interpretation of third-order statistics are possible 
and have high potential; our detection is in very good agreement with the WMAP7 best fit-model and this together with the fact that our 
measurement is robust against systematics, is a very encouraging outcome and an important proof of concept for future weak lensing missions. 
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